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1 Abstract 

Multi-stage time-stepping schemes that are tailored to chosen spatial-differencing operators 
are derived and tested. The schemes are constructed to give optimal damping of the high- 
frequency waves, making them ideal for use with multi-grid acceleration. The concept of 
characteristic time-stepping, necessary for the extension of the scalar analysis to systems 
of equations, is presented. The schemes show a marked improvement over Runge-Kutta 
schemes. 


2 Introduction 

Numerical methods for the computation of steady flows can be divided into two classes, 
explicit and implicit methods. Implicit methods have been favored for steady-state calcula- 
tions for a long time; this relates to the elliptic nature of the equations of steady subsonic 
flow. Representative of the class of implicit methods used for solving the steady Euler and 
Navier-Stokes equations are the Approximate-Factorization methods developed at NASA 
Ames Research Center [1] during the seventies. The strongest argument in favor of implicit 
methods is the relatively large reduction in residual that can be achieved in one iteration 
step. This is due in part to the numerical coupling of computational cells from one boundary 
to another in a single iteration. 

Explicit techniques require many more iterations than implicit techniques, but each itera- 
tion is relatively cheap. The numerical coupling from boundary to boundary can be achieved 
at low computational cost by the use of coarse-grid information, in a so-called multi-grid 
strategy. Once multi-grid relaxation has been successfully implemented, explicit methods 
have only advantages over implicit methods. They require little storage, are easily imple- 
mented on vector and parallel architectures, and naturally allow local grid refinements [2j. 
The latter advantage is crucial, since adaptive grid refinement seems to be the most promis- 
ing way to efficiently obtain spatial accuracy in complex problems. The most popular explicit 


1 



methods for computing steady solutions of the Euler and Navier-Stokes equations are the 
multi-stage methods pioneered by Jameson et al [3]. 

This lecture concerns the design of explicit multi-stage schemes for the Euler equations, 
for use in a multi-grid strategy. The requirements for single-grid and multi-grid schemes 
differ somewhat, although probably not as much as has been traditionally assumed. For 
single-grid computations the standard approach is always to take the maximum time-step 
allowed by the scheme s stability condition, the underlying idea being that the asymptotic 
steady state will then be reached in fewer steps. This is not necessarily the best strategy for 
multi-grid computations. 

In a multi-grid procedure, one special task of the marching scheme is to remove high- 
frequency components of the error while marching; the multi-grid strategy acts to remove 
low-frequency components through the use of coarse-grid representations of the solution [4], 
In all marching schemes presently in use, the best damping properties are achieved for a 
time-step that is substantially less than the maximum allowed by stability considerations. 

To make a marching scheme a good multi-grid “smoother,” the temporal and spatial 
discretization must be matched to each other. Since the spatial discretization dictates the 
final accuracy of the solution, the most natural way is to select a spatial discretization and 
design the time discretization in such a way that short waves are effectively damped. In 
the multi- grid code developed by Jameson [5], this is achieved by appropriately choosing the 
values of the parameters of a multi-stage integration method. The analysis on which the 
choice of parameters is based is strictly scalar and one-dimensional, and is carried out by 
trial and error. The material presented in this lecture suggests that a more thorough and 
comprehensive analysis of the damping properties of multi-stage schemes, and its proper ex- 
tension to systems of multi-dimensional equations, can significantly improve the performance 
of the multi-grid procedure. 

The core of this presentation is the analysis and optimization of the damping properties 
of multi-stage schemes for the one-dimensional linear convection equation (Section 2); the 
results can be applied to a nonlinear convection equation (Section 3), owing to the well-known 
scalar preconditioning technique of using “local” time steps. The extension to the system of 
the one-dimensonal Euler equations requires the use of “characteristic” time-steps, equivalent 
to preconditioning by a matrix (Section 4). The successful application to multi-dimensional 
scalar equations depends on the availability of a technique to damp numerical signals that 
move normal to the physical transport directions; application to the multidimensional Euler 
equations in addition requires a new matrix preconditioning. These two techniques are still 
in development (Section 5). 

The analysis is illustrated with numerical experiments throughout Sections 3, 4 and 5 
In Section 6 the results are summarized and a prognosis is given for explicit multi-grid 
relaxation for the multi-dimensional Euler equations. 
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3 The Multi-Stage Scheme as a High-Frequency Fil- 
ter 


The generic marching scheme used is a two-stage or predictor- corrector integration method 
for the linear ordinary differential equation 
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Here a is the time-step ratio and is a free parameter. As seen from Equation 2c, the stability 
and damping properties of the scheme are associated with the complex polynomial 


P -2 (z, a) = 1 + z + az 2 , z = AA< . 


(3a) 


This polynomial has two complex- conjugate roots, z\ (or) and Z 2 (a) — z \ ( a )> with 


(3b) 


these may be moved along the circle 
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(3c) 


by varying a. 

When a partial differential equation is interpreted by the method of lines, A represents the 
Fourier transform of the spatial differencing operator, and depends on the spatial frequency 
£ or, more specifically, on the spatial wave number 

/? = 2n(,Ax . (4) 


For instance, when solving the convection equation 

du du _ 

m = - c Tx' c>0 ’ 

use of upwind differencing for the spatial derivative gives 

1 

At— = — v [u (®, t) — u ( x — Ax, t )] , 


where the non-dimensional time step, 


v = 


cAt 
~Ax ’ 


(5) 

( 6 ) 
(7) 
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is the Courant- Friedrichs- Lewy (CFL) number. After inserting harmonic data 

ti(*) = u oe 2 ***, 

Equation 6 reduces to Equation 1 with 

XAt = z (0, v ) = —v (l — e - '^) . (9) 

The key observation to be made here is that, for any 0o in the high-frequency range 
[x/2, 7r], it is possible to make z(0o, v) = X(0o)At coincide with a zero ofP 2 (z, a), by choosing 
a particular combination of a and At. This results in perfect damping of the wave with wave 
number 0 O in one application of the predictor-corrector scheme. Using strings of predictor- 
corrector schemes, tuned to damp different frequencies, the entire high-frequency range can 
be damped to arbitrarily low levels. 

Strings of predictor-corrector methods generate multi-stage methods with an even number 
of stages; to get an odd number of stages, a single application of the “forward- Euler” scheme 

« n+1 = u n + AtXu n (10a) 

= (1 + XAt) u" , (10b) 

should be included in the string. The forward- Euler step has amplification factor 

Pi(z) = \ + z\ (11) 

this polynomial has one zero, at 

*i = -1 - (12) 

A second key observation is that, for any fixed number of stages, there is an optimum 
scheme, in the sense, reducing all of the high frequencies to an amplitude not exceeding a 
unique minimal threshold level. It is these “optimally smoothing” multi-stage schemes that 
axe developed here, for use in multi-grid Euler codes. 


3.1 Optimal Multi-Stage Schemes for General Spatial Differenc- 
ing 

For a general spatial-differencing operator, whether convective or diffusive, or a combination 
of both, the Fourier transform can still be written as 

XAt = z(0, i/) = v[a{0) + ib(0)\, (13) 

where u is a nondimensional time-step; keeping to the framework of convective equations for 
the sake of example, v shall continue to be referred to as the CFL number. For a frequency 
0o to be perfectly damped by the two-stage scheme, associated with the polynomial 3a, it 
is necessary to set z\{a) = z(0o, i/), or 

— + ;(<$>))] . (14) 
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with solution 



°o + bl 
°° “ 4 al ’ 

(15a) 


2 1« 0 1 

" '° - 

(15b) 

where a 0 = a((3 0 ), bo = Hj3 0 ). The single-step forward-Euler scheme 
quency for which the Fourier transform (Equation 13) is real- valued; 

can only damp the fre- 
for any finite-difference 

operator this means 

II 

(16) 

the corresponding CFL number is 

-U 

II 

(17) 


For the first-order upwind differencing operator, the optimization of the high-frequency 
damping in a string of predictor-corrector and single-step operators can be done analyti- 
cally [6]; for more complex differencing operators this is no longer feasible. An iterative 
method for solving this minimax problem is described below. 

Suppose that, in optimizing an m— stage method (m even or odd), as an initial guess or 
from the previous iteration, a set of perfectly damped frequencies /3 fc , k = l,...,m, with 
m = int(m/2), has been obtained; let these be represented by the vector (3 0 .^ There are 
m + 1 local maxima of |P m (/?)| on the interval [tt/2, tt] , called M k , k = 1, . . . ,m + 1, these 
can be found by a numerical search. As long as these values are not equal, the scheme is not 
optimal. The lack of equality is expressed in the form of an L 2 -residual 


m 

Po E: R{/3 0 ) = jT(M fc+1 (/3 0 ) - M k {(3 0 ))\ (18) 

Ar=l 

which is then brought closer to zero by one step of a Newton process. Since m frequencies 
must be updated, m independent residual values must be obtained. For this purpose, all 
frequencies f3 k are perturbed by a small amount <5/?, i.e. , 

& = & + 6/?, ( 19 ) 

and m new frequency vectors (3^ are formed such that the first i frequencies are perturbed, 

/3 W = ••>&.), * = (20) 

For all of these frequency vectors, residuals are obtained in the manner described above, 
these are called Pb) ? 

R.M = R(l 3 (,) ), i = l,...,m, (21) 

and are collectively indicated by the residual vector R : 

Hee(P ( 1 \...,p(-)). ( 22 ) 

The Jacobian matrix J = dR{p)/d(3 of the residual vector with respect to the frequency 
vector, needed for a Newton step, is computed approximately by means of finite-differencing. 
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Defining the vectors as f3^ with the _; th component either perturbed or ‘unperturbed’, 
i.e. 



= (ft , . . . A+l, • • .,£*), 

j < ft 

(23a) 


= (ftli • ■ • 5 fti /?•+!?• • • ifij-uftjyft+u • • • >/?m)> 

j > * * 

(23b) 

and defining the corresdponding residuals R^\ 




0$ 

III 

*5 


(24) 

some of which have been calculated previously; the elements of the Jacobian are 
uated according to 

then eval- 


_ dR (k) 

Ju = w, 


(25a) 


RW - /?<*•*> 

= v • l ^ k ' 


(25b) 


R{kA _ R (k) 

Jkt = xa , Z>k. 

bp 


(25c) 

The Newton step amounts to solving for the correction vector Af3 

from 



— JA/3 = R, 


(26) 

and updating (3 0 : 

(3 0 := (3 a +A(3. 


(27) 


The above procedure was implemented to find the first six multi-stage schemes with 
optimal high-frequency damping for the spatial differencing operator with Fourier transform 


\At = — i/(l — exp ,/? ) jl -| — — (1 — exp *^) H -j— (exp*^ — 1)| , (28) 

corresponding to higher-order upwind-biased differencing [7]. The parameter k regulates the 
upwind bias: k = 1 yields central differencing, k = — 1 second-order-accurate fully upwind 
differencing, k = 1/3 third-order-accurate upwind-biased differencing. For this latter choice, 
the damping properties of the six schemes are displayed in Figures 1 to 11. 

The important parameters of the first six multi-stage schemes for several spatial-differencing 
operators are listed in Tables 1-5; these require some explanation. Specifically, the quanti- 
ties ak listed in Tables 1-3 are not the time-step ratios of the constituent predictor-corrector 
schemes, but the time-step ratios arising in the practical implementation of an m-stage 


scheme, i.e., 

u<°> = u n , 


(29a) 


= u(°> T akAl , 

k = 1, . . . ,m , 

(29b) 


u n+1 = u (m) . 


(29c) 
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Magnitude of Amplification Factor for a = 1/3 



Figure 1: Amplification factor; third-order, one-stage scheme 


Magnitude of Amplification Factor for a = 1/3 



Figure 2: Amplification factor; third-order, two-stage scheme 
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Fourier Footprint of * = 1/3 
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Figure 3: Locus and contours; third-order, two-stage scheme 


Magnitude of Amplification Factor for * — 1/3 



Figure 4: Amplification factor; third-order, three-stage scheme 
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Figure 5: Locus and contours; third- order, three- stage scheme 


Magnitude of Amplification Factor for K — 1/3 



Figure 6: Amplification factor; third-order, four-stage scheme 
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Fourier Footprint of a = 1/3 



Figure 7: Locus and contours; third-order, four-stage scheme 


Magnitude of Amplification Factor for a = 1/3 



Figure 8: Amplification factor; third-order, five-stage scheme 
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Fourier Footprint of k = 1/3 
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Figure 9: Locus and contours; third-order, five-stage scheme 


Magnitude of Amplification Factor for k = 1/3 



Figure 10: Amplification factor; third-order, six-stage scheme 
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Fourier Footprint of k — 1/3 
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Figure 1 1 : Locus and contours; third-order, six-stage scheme 


Note that the m th step always spans the full time-interval At, so that 

«m = 1. (30) 

The amplification factor of the above scheme can be written as 

P m (z) = 1 + z(l + ct m -\z(\ + ot m -2z( (1 + a 2 z{l + ot\ z)) . . .))) , (31a) 

where z now corresponds to the CFL number for the full time- interval. The coefficients 
a* thus are found by multiplying out the string of polynomials of the form 3a and 11, and 
re-scaling z such that the linear term gets a coefficient of unity. 

For example, the three-stage first-order method, regarded as a string of a single step and 
a predictor-corrector operator, can be shown to have the amplification factor 


aw = (1 + §)(1 + * + § A 

(32) 

with £ corresponding to a CFL number of 1; it can be rewritten as 


J *>= 1+ ? + £‘ ,+ r’ 

(33) 

, /3 \ 2/3 \ 2 8 /3 \ 3 

+ \27 + 5 \2 / + 135 V2 7 

(33) 

or 


= 1 + * + !*’ + 155** 

(35a) 

= 1 +*('+!*(! + ^*)) . 

(35b) 
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Number of Stages 



1 

2 

3 

4 

5 

6 


1 

0.3333 

0.1481 

0.0833 

0.0533 

0.0370 



1 

0.4000 

0.2069 

0.1263 

0.0851 

0(3 



1 

0.4265 

0.2375 

0.1521 

<*4 




1 

0.4414 

0.2562 

« 5 





1 

0.4512 

0(6 






1 


Table 1: Multi-stage Coefficients for Optimal First-Order Scheme 


Number of Stages 



2 

3 

4 

5 

6 

<*i 

0.4242 

0.1918 

0.1084 

0.0695 

0.0482 

<x 2 

1 

0.4929 

0.2602 

0.1602 

0.1085 



1 

0.5052 

0.2898 

0.1885 

a 4 



1 

0.5060 

0.3050 

«s 




1 

0.5Q63 

a 6 





1 


Table 2: Multi-stage Coefficients for Optimal Second-Order (k = -1) Scheme 
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Number of Stages 



2 

3 

4 

5 

6 

<*1 

0.6612 

0.2884 

0.1666 

0.1067 

0.0742 

a 2 

1 

0.5010 

0.3027 

0.1979 

0.1393 



1 

0.5275 

0.3232 

0.2198 

a 4 



1 

0.5201 

0.3302 

a 5 




1 

0.5181 

<*6 





1 


Table 3: Multi-stage Coefficients for Optimal Third-Order ( k = 1/3) Scheme 


with z = (3/2)z. The scheme is completely defined by specifying 


4 

ai = 27’ 

(36a) 

2 

« 2 = 5, 

(36b) 

«3 = 1, 

(36c) 

3 

17 = 2- 

(36d) 


As mentioned before, the CFL number achieved in an m-stage scheme optimized for a partic- 
ular spatial-differencing operator is considerably lower than the maximum CFL number that 
can be realized using any m-stage scheme. For the first-order upwind-differencing operator, 
for instance, the maximum CFL number attainable in m steps equals m; in the scheme with 
optimized high-frequency damping the CFL number amounts to m/2. The minimax values 
of |P| in the high-frequency range are shown in Table 4; the CFL numbers to achieve these 
are shown in Table 5. 


4 Application to a Nonlinear Scalar Equation 

A nonlinear convection equation with a source term, 

t + JKir) = r in(,ri) ' ie[(MI ' (37) 

was chosen for a scalar test of the optimally smoothing multi-stage schemes. The spa- 
tial operator was approximated by third-order upwind-biased differencing, corresponding to 
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Number of Stages 



1 

2 

3 

4 

5 

6 

First 

Orrlpr 

0.7071 

0.3333 

0.1415 

0.0589 

0.0244 

0.0101 

«-i 


0.8093 

0.6521 

0.4309 

0.3030 

0.2073 

K = 1 3 


0.7016 

0.4668 

0.2950 

0.1848 

0.1153 

K = 0 


0.6636 

0.4213 

0.2579 

0.1558 

0.0940 

K = “3 


0.6432 

0.4009 

0.2413 

0.1435 

0.0851 

* = -§ 


0.6289 

0.3887 

0.2316 

0.1364 

0.0794 

AC = — 1 


0.6179 

. 

0.3801 

0.2244 

0.1315 

0.0759 


Table 4: \P\ max for 0 G [tt/2, tt] for Optimal Schemes 


Scheme 


Number of Stages 



l 

2 

3 

4 

5 

6 

First 

Orrlpr 

0.5000 

1.0000 

1.5000 

2.0000 

2.5000 

3.0000 

e^jro 

II 


0.9132 

2.0333 

2.3252 

3.0438 

3.5865 

*= 3 


0.8276 

1.3254 

1.7320 

2.1668 

2.5975 

K = 0 


0.7031 

1.0560 

1.3994 

1.7487 

2.0976 

K = ”3 


0.6055 

0.8950 

1.1885 

1.4844 

1.7802 

* = -§ 


0.5295 

0.7808 

1.0371 

1.2953 

1.5536 

K = — 1 


0.4693 

0.6936 

0.9214 

1.1508 

1.3805 


Table 5: Optimal CFL number for Optimal Schemes 
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Number of Stages 



2 

3 

4 

5 

6 

<*1 

1 

2 

I 

3 

1 

4 

1 

5 

1 

6 

a 2 

l 

1 

2 

1 

3 

1 

4 

1 

5 

<*3 


l 

1 

2 

1 

3 

1 

4 

a 4 



l 

1 

2 

1 

3 

<* 5 




1 

1 

2 






l 


Table 6: Multi-stage Coefficients for Runge-Kutta Schemes 


Equation 28 with k = 1/3. Steady solutions were sought on a grid of 512 cells. Two different 
kinds of marching-schemes were tried: 

1. Runge-Kutta multi-stage schemes, with use of the maximum stable CFL number; 

2. Optimally smoothing multi-stage schemes. 

When using the optimally smoothing multi-stage methods, it is crucial to make the CFL 
number constant over the entire grid, namely, equal to the unique value derived for maximum 
damping. This amounts to “local time-stepping” at the prescribed CFL number. For the 
Runge-Kutta schemes, local time-stepping was used at the highest stable CFL numbers. In 
cells where the convection speed passes through zero, the time-step must be limited (see, 
e.g. [8]). For both kinds of schemes, a saw-tooth cycle of multi-grid acceleration was used. 
All solutions were converged to a factor of 10“ 10 reduction in the residual. 

The coefficients for the first six Runge-Kutta schemes are shown in Table 6; the results of 
these schemes with regard to convergence speed are summarized in Table 7. The computa- 
tional work needed for convergence, expressed in terms of finest-grid residual calculations, is 
shown for various numbers of stages and grid levels. It is clear that the maximum-time-step 
strategy does not combine with the multi-grid strategy; the reason is that the Runge-Kutta 
schemes, like most schemes, are not good smoothers at the maximum stable CFL number. 

The results of the optimally smoothing schemes are shown in Table 8. These schemes 
clearly are a better match for multi-grid acceleration than the Runge-Kutta schemes. Con- 
vergence is reached more quickly in all cases investigated, even on a single grid, despite lower 
CFL numbers. It is interesting to note that the most efficient of the optimally smoothing 
schemes, for a sufficient number of grid levels, is the simple two-stage scheme, at least for 
this simple scalar problem. The gain in smoothing and CFL number achieved with a larger 
number of stages does not overcome the added computational work. 
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Number of Stages 


Grid 

Levels 



2 

3 

4 

5 

6 

1 

1352 

933 

816 

835 

942 

2 

921 

581 

528 

645 

585 

3 

739 

373 

350 

805 

641 

4 

1452 

315 

570 

919 

642 

5 

1585 

326 

721 

969 

744 

6 

2154 

378 

749 

1034 

839 


Table 7: Work required for convergence in scalar case — - Runge-Kutta schemes with maxi- 
mum stable CFL number. Work is expressed as number of finest-grid residual calculations. 


Grid 

Levels 


Number of Stages 



2 

3 

4 

5 

6 

1 

716 

627 

660 

655 

660 

2 

384 

365 

354 

360 

369 

3 

235 

221 

224 

237 

252 

4 

169 

180 

165 

188 

203 

5 

144 

198 

163 

175 

186 

6 

146 

201 

166 

178 

178 


Table 8: Work required for convergence in scalar case — Optimally smoothing schemes. 
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5 Application to a System of Equations 

The next series of numerical experiments was based on the quasi-one-dimensional Euler 
equations for flow in a converging- diverging channel. Only the optimally smoothing multi- 
stage schemes were tested. The equations were solved in the form 


(38) 





puA > 


< 0 \ 

d 

dt 

puA 

d_ 
' dx 

(pu 2 +p)A 

= 



\ pEA j 


K pu(E + p/p) A j 


\ 0 ) 


where A is the channel area, given by 


1 


A(x) = 1 + - (1 - COS (7TZ)) , 


-1 < X < 1 . (39) 

Two test cases were run, both with an inflow Mach number = 0.3059: 


1. Shockless transonic flow; 

2. Transonic flow with a shock in the diverging portion of the channel. 

The Mach number distribution for the second case is shown in Figure 12. In all cases, 
a sawtooth cycle was used for the multi-grid acceleration, with third-order upwind-biased 
spatial differencing (/c = 1/3) on the finest grid, and first-order upwind differencing on all 
coarser grids. In all cases, the parameters associated with the sawtooth cycle (number of 
solver applications on finest and on coarser grids) were varied in order to obtain convergence 
in the least amount of work. Each case was run with local time-stepping and characteristic 
time-stepping (explained below) for comparison. 

The choice of upwind-differencing is not coincidental; in fact it is mandatory for the 
analysis to extend to the system case. If a conventional approximation is used, based on 
central-differencing and an artificial viscosity with a scalar coefficient, the loci in the complex 
plane have a different shape for each wave mode. If upwind-differencing is used, the loci differ 
merely by a scale factor. Characteristic time-stepping, i.e. use of different time steps for the 
different characteristic equations such as to make all characteristic CFL numbers equal, can 
then be used to remove these scale factors. 

Characteristic time-stepping is equivalent to preconditioning the Euler residual by a 
local matrix, rather than a scalar, as in local time-stepping. To show this, the quasi-one- 
dimensional Euler equations are written as 


^ = - A (U)^ + s (U, x) = Res (U) ; 


(40) 


for the present analysis the conservation form is not required. The nominal matrix-preconditioned 
version of this equation reads 


^757 = max [p ( A )] I A l 1 Res > 

at x 


(41) 
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Mach Number Distribution in Channel 
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Figure 12: Mach number distribution — Transonic flow with shock 

where p( A) is the spectral radius of A, and |A| is the matrix with the same eigenvectors 
as A but with the absolute eigenvalues of the latter. In practice each of these eigenvalues 
may locally vanish, or be very small, making the inversion of |A| impossible or undesirable. 
Therefore, Equation 40 will actually be preconditioned according to 

^■7 = max[p(A)] (|A|*) _1 Res, (42) 

at 1 

where the eigenvalues of |A|* are bounded away from zero. The lower bound of the eigen- 
values of |A|* will depend on the magnitude of the local source term. The replacement of 
the matrix |A| by a non-singular matrix |A|* arises also in the numerical implementation of 
the so-called entropy condition for first-order hyperbolic system [8]. 

With only local time-stepping, the preconditioned Euler equations read 

^ = max [p (A)] [p (A)] -1 Res . (43) 

The work required for convergence of the shockiess cases, on a fine grid of 256 cells, 
is indicated in Tables 9 (local time-stepping) and 10 (characteristic time-stepping). For 
local time-stepping, one solver on the finest grid and one solver on the coarser grids was 
most efficient for all cases. For characteristic time-stepping, with six grid levels, it was 
advantageous to use two or three solver applications on the finest grid. Even then, the 
performance of the multi-grid relaxation leveled off for more than four grid levels. The 
reason for this is not yet understood. It is seen that characteristic time-stepping leads to 
a substantial improvement in convergence speed over local time-stepping. This is to be 
expected from the analysis, and can be traced to two causes: 

1. Characteristic time-stepping removes stiffness due to the variation among characteristic 
speeds, thus improving the performance of the scheme on a single grid; 
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Number of Stages 


Grid 
Levels 

Table 9: Work required for convergence — Shockless transonic flow — Local time-stepping 


2. By characteristic time-stepping, the optimal CFL number may be used for all waves 
simultaneously. 

The work required for convergence of the cases with a shock is shown in Tables 11 (local 
time-stepping) and 12 (characteristic time-stepping). For local time-stepping, the strategy 
of one solver application on each grid was the most efficient, except for a few cases with 
two or three grid levels. These cases required up to three solver applications on the finest 
and/or coarser grids. For characteristic time-stepping, the number of solver applications 
varied from one to eight on the finest grid, and from one to three on the coarser grids. Even 
with these variations, four of the cases did not achieve a residual drop of 10 -1 °, but got hung 
up between 10 -5 and 10 -9 . Nevertheless, in all converged cases, characteristic time-stepping 
helped speed up convergence. 

Figures 13-15 shed some light on the convergence properties of the schemes tested for 
the shock case. The comparison between local and characteristic time-stepping can be seen 
in Figure 13, which shows the two residual histories for the two-stage scheme with five grid 
levels. The effect of the number of stages on convergence can be seen in Figure 14, which 
shows the residual histories for characteristic time-stepping and five grid levels. The effect 
of the number of grid levels can be seen in Figure 15, which shows the residual histories for 
the two-stage scheme with characteristic time-stepping. 

That the multi-grid convergence is basically independent of the number of cells in the 
finest grid may be seen in Figures 16-19, which show the residual histories for a variety of 
cases. For each base grid, the numbers of stages and grid levels that gave the best perfor- 
mance were chosen. The convergence histories with characteristic time-stepping (Figures 17 
and 19) show the proper behavior, while for local time-stepping, shown in Figures 16 and 18, 
there is some dependence on the number of cells in the base grid. 
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Number of Stages 


Grid 

Levels 



Table 10: Work required for convergence — Shockless transonic flow Characteristic time- 
stepping 


Number of Stages 



Table 11: Work required for convergence — Transonic flow with shock — Local time-stepping 
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Number of Stages 


Grid 

Levels 



2 

3 

4 

5 

6 

1 

6130 

5200 

5107 

— 

4584 

2 

1589 

1608 

— 

— 

— 

3 

916 

863 

987 

850 

851 

4 

639 

536 

526 

587 

645 

5 

374 

365 

390 

332 

380 

6 

627 

429 

464 

482 

579 


Table 12: Work required for convergence — Transonic flow with shock 
time-stepping. Dashes denote cases that did not converge to lO -10 . 


Characteristic 


Convergence History 
Channel Flow 



Local 

Char. 


Figure 13: Comparison of local and characteristic time-stepping — Transonic flow with shock 
— Two stages, five grid levels 
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Convergence History 


Channel Flow 



2 Stages 

3 Stages 

4 Stages 

5 Stages 

6 Stages 


Figure 14: Convergence with different numbers of stages — Transonic flow with shock — 
Characteristic time-stepping, five grid levels 


Convergence History 
Channel Flow 



Figure 15: Convergence with different numbers of grid levels — Transonic flow with shock 
— Characteristic time-stepping, two stages 
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Convergence History 


Channel Flow 



Figure 16: Convergence with different base grids — Transonic flow without shock — Local 
time-stepping, best results per base grid 


Convergence History 
Channel Flow 



256 Celia 

512 Celia 

1024 Celia 


Figure 17: Convergence with different base grids — Transonic flow without shock — Char- 
acteristic time-stepping, best results per base grid 
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Convergence History 
Channel Flow 



64 Cells 

128 Cells 

.... 256 Cells 


Figure 18: Convergence with different base grids — Transonic flow with shock — Local 
time-stepping, best results per base grid 


Convergence History 
Channel Flow 



128 Cells 

256 Cells 

.... 512 CeUs 


Figure 19: Convergence with different base grids — Transonic flow with shock — Charac- 
teristic time-stepping, best results per base grid 
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6 Application to Multi-Dimensional Equations 

A successful extension of the scalar, one-dimensional analysis to the Euler equations in more 
than one space- dimension stands or falls with the availability of a robust wave-decomposition 
model. Work on such decompositions is in progress, but has not yet led to reliable schemes. 
The basis for one possible preconditioning follows. 

The two-dimensional Euler equations can be written as 


dV 

dt 


= _ A ( U)^- B( U)^=Res(U) 


(44) 


There is no obvious way to precondition the residual with a local matrix, as the matrices 
A(U) and B(U) do not have the same eigenvectors, and therefore can not be diagonalized 
simultaneously. This means Equation 44 can not be written as a system of coupled scalar 
convection equations. An understanding of the waves that locally pass through the grid can 
be obtained by assuming some a priori knowledge about the types of waves present, and 
then fitting this model to the local data [9, 10]. The wave model adopted here is based on 
the Euler equations written in a coordinate system aligned with the local streamline: 


dU 

dt 


= -A,(U)^-A„(U)^=Res(U). 


(45) 


The four-component residual vector and, in addition, the two components of the cell-averaged 
pressure gradient, are then used to obtain the amplitudes of six local waves: two acoustic 
waves, a shear wave and an entropy wave, all moving along the streamline, and two acoustic 
waves moving normal to the streamline. Mathematically this means that the residual is 
rewritten as 

Res (U) = RAa , (46) 

where a is the vector of wave strengths, A is a 6 x 6 diagonal matrix carrying the wave 
speeds in its main diagonal, and R is a 4 x 6 matrix built from four eigenvectors of A, and 
two of A n . The preconditioned equation may now formally be written as 


dU 

dt 


Kp (A) A -1 L Res (U) 
M Res (U) . 


(47) 

(48) 


The 6x4 matrix L is a generalized inverse of R and therefore not uniquely defined. Research 
is presently focusing on finding an inverse, using physical or mathematical arguments, that 
makes the preconditioning by M truly effective. 

It should be noted that, even if a suitable preconditioning matrix M is developed, the 
optimally-smoothing schemes are “tuned” to filter waves travelling in the propagation di- 
rection ( c x ,Cy ) of the problem being solved. For optimal performance in two dimensions, 
something must be done about high-frequency waves normal to this direction. For a convec- 
tion problem defined by 


d d 

Cx dx + Cy dy ’ 


(49) 
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Stages 


Grid 

Levels 


Table 13: Work required for convergence — Two-dimensional case — Local time-stepping 



2 

6 

1 

8821 

7404 

2 

2455 

3145 

3 

2458 

1758 

4 

2380 

1684 


a derivative in the direction normal to the propagation direction is given by 

d_ _ 3_ 

Cy c> a i 

ox ay 


(50) 


and a positive-definite operator 


which acts in this normal direction is given by 


2 d 2 u d 2 u 2 d 2 u 

Cx 7fh, ~ lCxC ’‘ai<h, + ' 


(51) 


Unfortunately, adding this term would reduce any scheme to first order. Methods of adding 
this term in a non-linear manner are being studied, to give optimal damping without reducing 
the order of the solution in the steady state. 

Since the work on the preconditioning matrix and the cross- diffusion terms mentioned 
above is still in progress, the two-dimensional Euler equation results presented here are for 
local time-stepping with no cross-diffusion. The test case was a NACA 0012 airfoil at zero 
incidence in a M <*, = 1.2 freestream. This case has a bow shock and a fishtail shock. First- 
order upwind-differencing was used. Table 13 summarizes the number of finest-grid residual 
calculations necessary for convergence (to 10 -1 °) on a 64 x 32 grid. In the six-stage, four 
grid-level scheme, three iterations on the finest grid followed by one on each of the coarser 
grids were used; for all other schemes, one iteration on each grid was used. Despite the use of 
local time-stepping as a substitute for characteristic time- stepping, the multi-stage schemes 
combined well with the multigrid acceleration. 


7 Conclusions and Future Research Directions 

In these notes, a method has been developed for designing optimally smoothing multi-stage 
time-marching schemes, given any spatial-differencing operator. Such schemes are particu- 
larly useful in conjunction with multi-grid acceleration. The advantage of using these opti- 
mally smoothing schemes has been demonstrated by comparison with Runge-Kutta schemes 
in solving a nonlinear scalar equation. The analysis has been extended to the Euler equations 
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in one space-dimension by use of characteristic time-stepping. Convergence rates indepen- 
dent of the number of cells in the finest grid have been achieved with these optimal schemes, 
for transonic flow with and without a shock. Besides characteristic time-stepping, local 
time-stepping has been tested with these schemes. While the analysis is only truly appli- 
cable with characteristic time-stepping, good convergence has still been obtained with local 
time-stepping. The extension to two-dimensional flows is hampered by the lack of a robust 
two-dimensional wave model that may serve as the basis of characteristic time-stepping, and 
by the lack of a method to damp high-frequency waves normal to the direction of propaga- 
tion. Future research must concentrate on these two issues. Only with these techniques may 
full advantage be taken of the optimally smoothing multi-stage schemes. 
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